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ABSTRACT 

The objective of this paper is to describe an accurate and 
efficient reduced order modeling method for aero elastic (AE) 
analysis and for determining the flutter boundary. Without los- 
ing accuracy ; we develop a reduced order model based on the 
Volterra series to achieve significant savings in computational 
cost. The aerodynamic force is provided by a high-fidelity so- 
lution from the Reynolds -averaged Navier-Stokes (RANS) equa- 
tions; the structural mode shapes are determined from the finite 
element analysis. The fluid- structure coupling is then modeled 
by the state-space formulation with the structural displacement 
as input and the aerodynamic force as output , which in turn acts 
as an external force to the aero elastic displacement equation for 
providing the structural deformation. NASA’s rotor 67 blade is 
used to study its aeroelastic characteristics under the designated 
operating condition. First, the CFD results are validated against 
measured data available for the steady state condition. Then, 
the accuracy of the developed reduced order model is compared 
with the full-order solutions. Finally the aeroelastic solutions 
of the blade are computed and a flutter boundary is identified, 
suggesting that the rotor, with the material property chosen for 
the study, is structurally stable at the operating condition, free of 
encountering flutter. 
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NOMENCLATURE 

a Speed of sound. 

A Area vector. 

c Blade chord length. 

C p Specific heat at constant pressure, 
d structural deformation vector. 

E Specific total energy. 

F m Modal force vector. 

F Inviscid flux vector. 

F v viscous flux vector. 

h m (n) m^- order Volterra series kernel. 

K Stiffness matrix. 

L Length measured in the x-direction. 

M Mass matrix. 

Moo Freestream Mach number. 
p Pressure. 

Pk(co) Turbulence production in the k(co) equation. 

Pr Prandtl number. 
qi Heat flux vector. 

qoo Dynamic pressure(=Pootfd) 

U Conservative variables. 
r Position vector. 

Re Reynolds number. 

S Source term. 

T Temperature. 

Uoo Inflow velocity magnitude. 

ii[ Cartesian velocity components. 

x,y,z Axial (horizontal), span wise and vertical direction. 

Sij Kronecker delta. 

£ , rj Modal displacement. 
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K Turbulence kinetic energy. 

/I Viscosity. 

CO Specific turbulence dissipation rate. 

Q Rotating speed. 

a Value of step function. 

T Pseudo time. 

T ij Viscous stress tensor. 

"V Control volume. 

Subscripts 

A Aerodynamic. 

°o Far upstream or “free” stream. 

m Modal coordinate 

s Structural dynamics. 

t Total (stagnation) condition. 
turb Turbulence. 


1 INTRODUCTION 

NASA is considering new generations of aircraft that meet 
aggressive economic, noise and environmental targets; a spe- 
cial configuration, called N3-X, employs all electric power and 
propulsion systems by which the thrust force for the vehicle is 
generated exclusively with an array of fans housed in compart- 
mentalized flow paths. Hence, the designing of fans to meet es- 
sential considerations is paramount. Specifically, increasing per- 
formance and operating life and reducing weight to optimize the 
economic objective, while reducing noise and emission to meet 
environmental regulations. In pursuit of higher performance of a 
compressor/fan, the past design trend is to run at higher pressure 
ratios and higher mass flow rates, thus moving close to flutter 
boundaries associated with surge or choke as defined in the com- 
pressor map. Hence, it is important to ensure the compressor 
is structurally sound over the entire operating range, from the 
choke to the stall conditions. Structural vibration, either caused 
by natural resonance or forced response, is a major consideration 
in assessing the devise’s structural integrity. The fluid-induced 
instability of a compressor blade is typically not of concern un- 
less it is tuned to the natural vibration frequency. However, it 
becomes an issue in transonic speed regime, because a small dis- 
turbance can result in a large amplitude variation and nonlinear 
behavior. The unsteady excursion of a shock wave through the 
blade-to-blade passage can intermittently choke or stall the flow, 
potentially crossing the flutter boundary. The unsteady forces 
resulting from the shock motion are shown to have either stabi- 
lizing or destabilizing effects, depending on the shock structure 
and inter-blade phase angle. [1] 

With advances in computers and computational fluid dynam- 
ics (CFD), aeroelastic analysis is fast becoming common for real 
world designs. To be useful and adopted in practice, a compu- 
tational tool must be reliable for predicting the aeroelastic char- 
acteristics and just as importantly be efficient (cheap and fast). 
This computational tool will consist of an aerodynamics code 


and a structural dynamics code, as a result they in turn determine 
the reliability and efficiency of the tool. 

For aerodynamics analysis, developments in computational 
fluid dynamics over the past several decades have provided in- 
creasingly powerful and reliable capabilities. The complex- 
ity and fidelity, hence its range of applicability, of analysis is 
strongly correlated with the fast evolution of computer power: 
from the early linearized potential flow solution to the current 
large eddy simulations using Navier- Stokes equations. Linear 
models are still used widely in the design phase. But develop- 
ments in computer technology and CFD methods and software 
have made use of high-fidelity models feasible even in early stage 
of a design cycle. However, large eddy simulations are still far 
too costly and from being timely to be adopted in the design pro- 
cess. 

In the current study, we employ the Reynolds-average 
Navier-Stokes equations for which the turbulence terms are 
closed with the two-equation K-co model, specifically the shear 
stress transport (SST) version by Menter [2]. The second-order 
backward differencing is used for time-discretization. The non- 
linear inviscid terms are approximated by the AUSM + -up [3] 
method while the viscous terms approximated by the usual cen- 
tered formulas. The resulting implicit algebraic system is then 
solved by the LUSGS method [4]. 

For structural dynamics, one may invoke the full finite ele- 
ment analysis, as employed in the aeroelastic study of rotor 67 
by Doi [5]. The resulting fluid- structure system is a time depen- 
dent set of equations describing not only the flow variables in the 
entire domain, but also the motion of the structure immersed in 
the fluid. The system can be solved either in the frequency [6] or 
time domain [5]. The frequency domain approach may be pre- 
ferred for linear problems for its computational efficiency; how- 
ever for a nonlinear problem, it is more efficient and accurate to 
arrive at solution with the time domain approach. 

The time-domain computation for flutter analysis can be- 
come costly when a large number of time-dependent solutions 
of the fluid-structure system are needed, for example as part of 
a design process. It is therefore desirable to reduce the com- 
putational cost by a significant factor, for example by at least 
an order of magnitude or more. This can be readily achieved 
by employing strategies called model order reduction (MOR), of 
which the harmonic balance, proper orthogonal decomposition, 
and Volterra series are among the most popular. [7] Model order 
reduction should not only save computational effort, but also re- 
tain the fidelity of the original (full) system. This goal has been 
well realized for linear problems through model reduction, but 
not yet universally for nonlinear problems. 

For nonlinear problems, the Volterra series expansion is used 
to approximate the input-output relationship of a nonlinear time 
dependent system, with a capability of capturing ’’memory” ef- 
fects. This input-output concept is well suited for the aeroelastic 
analysis in which the aerodynamic force and structural deforma- 
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tion can be formulated in this framework. Furthermore, flowfield 
and structural dynamics have different time scales and their in- 
teractions often respond with delay in time, i.e., with ’’memory” 
effects. The Volterra series has been applied in various fields 
of engineering and is mostly used to construct a reduced order 
model to mimic a complex dynamic system. Unsteady aerody- 
namic force responses to wing motion have been calculated by 
Silva [8] using the Volterra theory. In the present work, we de- 
scribe the application of the Volterra series, based on RANS so- 
lutions, to turbomachinery aeroelastic problems. 

The paper is organized as follows. Section 2 gives the 
fluid and structure equations employed and outlines the methods 
adopted to solve them, especially including detailed description 
of fluid- structure coupling and model order reduction based on 
the Volterra series. Section 3 presents the application to aeroelas- 
tic analysis for NASA’s rotor 67 compressor blade along with the 
validation of the CFD solution against measured data. In Section 
3.2 we show the results of applying the developed reduced-order 
model to find the flutter boundary of rotor 67. 
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2 MATHEMATICAL SYSTEM 

The mathematical system considered for predicting the flut- 
ter conditions consists of the fluid and structure equations, as 
described in two respective sections. 


Here the subscript “/” denotes the direction in Cartesian coordi- 
nates, e.g., F = Y^=\ F/ e /. For turbomachinery applications, the 
coordinate x\ is chosen to be aligned with the rotor axial direc- 
tion, and the grid velocity follows a rigid body rotation with Q/, 

u g = rx£li = Q(zj- yk) (4) 


2.1 Fluid Equations 

The three-dimensional Reynolds-averaged Navier-Stokes 
(RANS) equations are employed, with the turbulence described 
by the two-equation zc-O) SST model [2]. They are written in the 
following integral form over a control volume Y(x,t) enclosed 
by a control surface dY (x,f): 


d 

dt 


f i 

Jy 


U dV- 



(F + P) • dA 


l F v • dA + [ S dV 
Jdy Jy 


(i) 


where we have the standard notation for the conservative vari- 
ables plus the turbulence variables in U. The surface integral on 
dY (x,t) consists of fluxes through the vectorial area dA can be 
expressed in terms of 3 Cartesian coordinates. The relative con- 
vective flux F i, the pressure flux P;, and the viscous stresses and 
heat flux FV in the /-direction, / = 1,2,3 are given in Eq. (2-3), 
written in the relative coordinate system moving with the speed 
u g . [9] The source terms includes the rigid-body rotation, Eq. (4) 


It is noted that in the aeroelastic calculation performed in this 
study, the position vector r(t) of a computational cell varies with 
time as the blade geometry changes, even though the rotating 
speed remains fixed. 

The above viscous stress terms include both the laminar and 
turbulent effects through the use of eddy viscosity (Ht,k t ) model, 
expressed as: 


% = (P + AW(>) 

(E I Bur£ ) 3 t_ 
qi \ Pr ^ Pr turb ) dxi 


duj | duj 2 du k x 
dxi + dx x 3 dx k °' J 


( 5 ) 


The turbulence eddy viscosity for the study presented herein 
is provided through the solution of two transport equations of 
scalar quantities (zc, co), specifically the so-called k-co turbu- 
lence model [10] is enhanced with Menter’s shear-stress trans- 
port (SST) model [2]. The details of the turbulence model, well- 
known and elaborated in the cited reference, are omitted here. 

2.2 CFD Solution Methods 

To eliminate accumulative time integration error, we opt for 
the dual-time stepping approach, in which a time rate of change 
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of U in pseudo time (t) is added to Eq. (1) so that at each 
new time level the unsteady CFD equations are balanced, namely 
driving its discretized residual to diminish. The dual time step- 
ping strategy is formulated as: 


— / U dV = R(U) 
dx Jy 

= ---[udV-(£ (F + P) • dA + (f F v -dA+ [ SdX6) 
dt Jy Jdy Jdy Jy 

In the physical time step (t) the residual R(U) is discretized 
implicitly using a 3-level, backward differencing in order to ob- 
tain second-order temporal accuracy, resulting in a highly nonlin- 
ear system. In the pseudo-time step the implicit system is solved 
by performing fixed-point iterations untill the residual R(U) of 
the nonlinear physical-time equation has diminished or reduced 
to specified small values. Then the solution is advanced to the 
next time level. This pseudo-time iteration is carried out by em- 
ploying the LUSGS method [4]. 

The inviscid flux terms F + P arguably have received the 
most attention in past CFD research, especially those under the 
framework of upwind solvers, yielding many proposed schemes 
for approximating it. In this study, We employ the AUSM + -up 
method [1 1] for the inviscid fluxes. For the viscous terms F v and 
source S terms, a typical centered representation is used. 

The mesh velocity is obtained from the structural motion in 
response to the aerodynamic forces provided by the CFD solu- 
tion. The structural model is described in the next section. 

The resulting in-house 3D RANS code has been developed 
and validated for a variety of flow problems over a number of 
years. For the validation relevant to the problem at hand will be 
described Section 3.1. 

2.3 Structural Dynamics Equation 

The finite element model for describing a structural motion 
is expressed in terms of its displacement t; from a neutral posi- 
tion (steady state in our case). In our work, we first carry out 
finite element analysis on a given set of nodes via MSC/Nastran 
[12] to obtain mode shapes, <£>/, i = 1,2 , • - - ,2V m , N m being the 
number of modes. Neglecting damping, the structural motion in 
terms of the modal displacement vector § = = 1,* • • ,N m } 

in response to the modal force F m can be described by 

M| + K§=F m (7) 

where (M,K) are the mass and stiffness matrices of the mate- 
rial of the structure respectively. The modes on the FEM nodes 
are then interpolated to every CFD node at which the aerody- 
namic forces are known. This modal information facilitates the 


determination of the modal (generalized) force F m = {F mi ; i = 
1, • • • ,N m }, where each component of the modal force is the in- 
ner product of the mode shape and aerodynamic force vectors 
over the entire CFD nodes. Additionally, the physical deforma- 
tion d of a structure can be expressed in terms of the mode shape 
and the modal displacement £ from the above dynamics equa- 
tion, leading to physical displacement, d = and physi- 

cal force f = <Z>,F mr 

The above second-order differential equation can be recast 
into the following first-order differential system: 

V = AjJ] +B s F m (8) 

where 
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and 77 = {?],■;/ = 1 

The time derivative in Eq. (8) is approximated by the 
second-order Crank-Nicolson method, producing a discrete sys- 
tem for t n <t < t n+x , 

V ( n + 1) = (I — 0.5Af A,)” 1 ((I + 0.5 AfA s )t] (n) + B sd F m (n )) 

( 10 ) 

This will form part of the coupled fluid- structure (aeroelas- 
tic) system to be elaborated below. It specifically provides the 
time-dependent modal displacement £ , hence the needed phys- 
ical displacement of the structure so as to affect fluid flow in 
response to the geometry variations. The mathematical system 
describing the interactions between fluid and structural dynam- 
ics is given below. 

2.4 Fluid-Structural Coupling 

The coupling of aerodynamic and structural computations 
must be performed on a common geometry, while they need not 
be of the same mesh density or matching at the same grid points, 
as displayed in Fig. 1 for the NASA rotor 67 blade, which is the 
structure that will be considered in this paper. As such, interpo- 
lation/extrapolation procedures must be employed to accomplish 
the mapping between them, through which the proper transfer of 
relevant variables may be carried out. In our case, the structure 
deformation provides a new body to the CFD process, thus af- 
fecting boundary condition and the flow domain mesh. On the 
other hand, the aerodynamic force needs to be transferred to the 
contact points for the finite element analysis. This mapping of 
grids must satisfy certain physical requirements, such as conser- 
vation of virtual work, and numerical requirements of accuracy 
and stability. 
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FIGURE 1: ILLUSTRATION OF GRIDS USED FOR 

AERODYNAMIC AND STRUCTURAL ANALYSES. BLUE: 
STRUCTURAL GRID; RED: AERODYNAMIC GRID. TYPI- 
CALLY THE STRUCTURAL GRID IS MUCH COARSER. 


The constant volume transformation method [13] was at- 
tempted to interpolate/extrapolate between the aerodynamic and 
structural grids, but it failed to provide a stable and converged 
solution because of severe geometrical twisting involved in the 
present case. Instead, we employed a surrogate model to provide 
the structural mode shapes. 

The radial-based function (RBF) neural network method, 
which we previously used for reduced-order modeling of flut- 
ter and limit cycle oscillations [14], is applied here by taking the 
finite element nodes as input and modes as output via a training 
process by virtue of Nastran calculations. Once the training (the 
RBF and neurons) is accomplished, the neural network can take 
the aerodynamic grid as input and produce modes as output, the 
modes are in turn interpolated to the aerodynamic grid. 

The resulting first three mode shapes are displayed in Fig. 
11, revealing indeed a large deformation of the structure. Once 
the blade shape deformation caused by the aerodynamic force 
is updated as described above, then the computational mesh for 
CFD is changed accordingly using transfinite interpolation (TFI) 
at each time step. As such, it also allows physical variables to be 
interpolated onto the new grid in the same manner. 

A typical fluid-structure coupling is performed as shown in 
Fig. 2, where the structural deformation is known, hence mesh 
generated at t n , a subsequent CFD solution for the new time step 
at t n+l for U" +1 is performed with the structural shape frozen 
at q n , then it is followed by a geometry update to get q n+1 
with the input of U^ +1 by procedure 3. This is the so-called 
loosely coupling strategy, in contrast to the tightly coupling one 



FIGURE 2: CFD AND CSD INTEGRATED COMPUTATION: 
LOOSE COUPLING. 


in which both the CFD and computational structural dynamics 
(CSD) equations are solved simultaneously. The loosely cou- 
pling strategy is easy to implement and computationally effi- 
cient; the lag in time between CFD and CSD is believed to be 
insignificant, since the time step is usually much smaller than the 
characteristic time of the problem under study. This combined 
CFD-CSD full order modeling process is extremely expensive 
especially when a large number of computations are committed. 

In real-world engineering practice, analysis is not performed 
only for one condition, but over many computations. In ad- 
dition, design optimization typically will require hundreds and 
thousands similar computations, differing for example in range 
of conditions, parameters, or geometry. A reduced-order model 
takes only a small fraction of computational time needed by a 
full-order model, but is of value if and only if it is capable of 
preserving the accuracy of the full-order system. This can be 
achieved easily for a linear system, but still remains a topic of 
intensive research for a nonlinear system [7]. 

As all the aeroelastic computations for finding flutter/LCOs 
are similar in kind and repetitive, they differ only by a limited 
number of variables and the variations in value. A model order 
reduction will be of great value in significantly reducing com- 
putational cost and time. In what follows we will describe the 
application of Voterra theory [15] for constructing a reduced or- 
der model for aeroelastic analysis, based on the fidelity of the 
RANS equations for aerodynamic calculations. 

2.5 MODEL ORDER REDUCTION BY VOLTERRA SE- 
RIES 

The Volterra theory provides a functional relationship repre- 
senting a nonlinear response to a given input function which may 
be time dependent and is capable of capturing the ’’memory” ef- 
fect. While it has been employed in previous studies for aeroe- 
lastic application, for example [16, 17], these have been limited 
to external flows over an airfoil or a wing. To our knowledge, the 
current paper represents the first aeroelastic application of the 
Volterra theory to turbomachine. The Volterra theory has some 
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advantages over other ROMs, see [7, 17] for more discussion. 
The Volterra theory can be easily adopted as an alternative proce- 
dure without having to modify the baseline full-order procedure. 
It is equally applicable to the time and frequency domains and 
the conversion between them is rather simple. Moreover, the for- 
mulation facilitates to retain nonlinearity of the full order model 
more easily than other ROMs. 

The Volterra series, unlike the Taylor series, includes in the 
output accumulative effects of inputs occurring at previous times. 
The output y(t) of a continuous time-invariant system in response 
to a single input u(t) for t > 0 is expressed by the Volterra theory 
as: 


a system response after applying a step function, 

(i4 > 

A small number is given to Go = 1.0 x 10 -4 to ensure the prob- 
lem remains linear. Then, according to Eq. (13), we have the 
response, 

n 

y(n) =h 0 + o 0 J^h \(n-k), (15) 

k=0 


y(t) = h o + 


if-fut 

£ 1-/0 Jo 


- Ti ,•••>*- 1 ) • • • u(ti)dti ■ ■ ■ dti. 


(ID 

where ho is the steady- state term coincident with the initial con- 
dition and hi , i > 1 are known as the Volterra kernels. As the time 
integral is discretized over a n-interval domain, a time-discrete 
infinite (or truncated) Volterra series is obtained: 


n 

y{n) = ho + ^ h\ (n — k)u(k) 
k = o 

n n 

+ Y Y h2(n — k\,n — k2)u(ki)u(k2) ( 12 ) 

= 0 * 2=0 
n n 

+ £•••£ h m {n-ki,---,n — k m )u(k\ )■■■ u(k m ) H 

k\ =0 k m = 0 


where y(ri) is the output with the time index n referring to t n , 
u(k) is the input at preceding times k = 0, 1,2, ...,n, and h m the 
rath-order Volterra kernel, ra = 1 , 2, • • • , ©o. For a linear system, it 
suffices to keep only the first-order Volterra kernel, hence 


n 


y(ri) = ho+ ^h\(n — k)u(k) 
k = o 


(13) 


where ho corresponds to the response with zero input, or the force 
vector at steady state where there is no structural response. To 
capture behaviors varying with time variation, one must at least 
find the first-order kernel associated with the input at all other 
times, it turns out that from the continuous system, the first ker- 
nel measures the response to an impulse applied at T\ = 0. To 
include nonlinear effects, higher order kernels are necessary, see 
Silva [8] 

In the present study, we make use of the first kernel to build 
our reduced order model (ROM), for which the necessary step is 
the definition of h\ (n), for n > 0, as will be illustrated below for 


And the first kernel is readily available as 

hl{n)= { 0 (kn)-y(n- l))/oo, «>L (16) 

The first equality holds because of the initial condition y(0) = 

K o). 

In what follows we show how to construct a reduced order 
model that simply bases on a relationship between the structural 
motion and aerodynamic force, from the viewpoint of relating 
input and output data. This is easily facilitated within the state- 
space theory, as used in control theory, a linear state-space sys- 
tem can be represented in the following canonical form: 

Xfl(n-bl) = A fl x fl (ft)+B^(ft) (17) 

F a (rc+ 1) — C a x a (n) -\-D a t;(n) (18) 


where x a (n) is the state vector at time n. The input t, is the 
structural displacement and the system output ¥ a denotes the non 
dimensional generalized aerodynamic force. 

To set up the above system and solve for the aeroelastic sys- 
tem under consideration, we adopt the Eigensystem Realization 
Algorithm (ERA) [18]. First, we define the finite Hankel matrix 
constructed using the first-order Volterra kernel Ai just described 
above, 


H(k- 1) 


h i(fc) 

hi(k+ 1) 

h\ (k-\- p — 1) 

hi(k+ 1) 

hi(k + 2) 

hi(k + [i) 

hi(k + 2) 

hi(k + 3) 

h\ (k + j6 + 1) 

h\ (k ~t - oc — 1) 

h\ ( k -f Of) • 1 

• • h\ {k oc -f J3 — 


where a and ft are the sampling time shift in the row and column 
directions respectively; they control the order (rank) of the sys- 
tem, and are set as a = 1600, and /3 = 50 in our study. Applying 
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Singular Value Decomposition to H( 0), 

H( 0) = ULV t (20) 

we find £/, E and V, which are then used to define the matrices in 
Eq. (17): 

A a = E - 1/2 U T H(0)VL~ l/2 
B a = I}! 2 V t E l 
C a = E T M Ul}! 2 

D a = h\ (0) (21) 

where 

e m = [Im 0m ■■■ 0 M ] aMxM 

E t l =[Il 0l ■■■ 0l] PL xL ( 22 ) 

with M and L being the number of inputs and outputs respec- 
tively. Since only first three modes are retained, we have M = 3 
and L = 3. The size of the ROM is 3 x j 3 = 150 for this study. 



n 

FIGURE 3: DECAY OF (NORMALIZED) SINGULAR VAL- 
UES OF THE HANKEL MATRIX. 


In Fig. 3, we show the efficiency of the reduction method. 
Singular values of the Hankel function are seen to decay rapidly 
within the first 50 values, indicating that model order is of rea- 
sonable size needed to retain accuracy. 

It is appropriate at this juncture to illustrate the entire aeroe- 
lastic analysis process in a flow chart shown in Fig. 4. The ROM 


track starts with the baseline CFD solution as the full-order will, 
then builds the Volterra kernel shown in Eq. (16), which forms 
the state space system in Eq. (17). The input and output of which, 
t, and F a , are coupled with the structural dynamics system in Eq. 
(10). It is noted that F m = with being the dynamic pres- 
sure ( poodl , ). These two systems combined form the ROM for the 
aeroelastic analysis discussed next. 



FIGURE 4: FLOW CHART ILLUSTRATING THE PRO- 

CESS OF PERFORMING AEROELASTIC ANALYSIS IN 
THE PRESENT STUDY. 


3 Aeroelastic Analysis of NASA Rotor 67 

Complex vibration problems arise from the interactions of 
nonlinear aerodynamics and structural deformation. These vi- 
brations can be either self-induced or caused by flow distortions 
from upstream and downstream blade rows or tip region; the for- 
mer is called flutter and the later forced response. Also mistuning 
in blade rows and inter-blade phase angle can force vibration on 
a blade. [19] Bendiksen [20] gives a comprehensive review of 
aeroelastic problems encountered in turbo machines, in which 
various factors causing flutter are identified. A recent discus- 
sion on the progress and challenges of computational aeroelsac- 
tic modeling can be found in Bartels and Sayma. [21] Aeroelas- 
tic analyses of rotor 67 have been conducted by Doi [5], Sadeghi 
and Liu [22] and Zhang et al [19], employing full-order mod- 
eling of the fluid- structural system where RANS is used in the 
CFD procedure. An inlet guid vane row is included in Zhang 
et al to study its effect on the flutter characteristics. Doi found 
that operating condition and inter-blade phase angle determine 
the stability of the structural response; Sadeghi found the rotor 
to be stable using the 10 first eigenmodes. In our study, we also 
choose this rotor for our computational test model because the 
model reveals most flow complexities seen in the turbomachines 
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in today’s aircraft and also widely used in the turbo machinery 
community for validation of CFD results, thus allowing us to 
verify our proposed approach for AE analysis against previous 
works, for example [23, 24, 5, 22]. 

In what follows, we shall first validate the CFD solution for 
detailed profiles and performance map against the measured data 
taken in [25]. Then the fluid- structure coupling procedure will be 
described, followed by the aeroelastic calculation of the blade. A 
model order reduction method based on the Volterra series is in- 
troduced and applied to rotor 67 to determine the flutter behavior. 

3.1 Validation at Steady State Operating Points 

NASA rotor 67, shown in Fig. 5 is the first stage rotor of a 
two-stage fan [26]. It is a low aspect ratio (1.56) transonic axial 
flow rotor with a design tip relative Mach number of 1.38; an ex- 
periment program was undertaken to provide laser anemometry 
and aerodynamic performance data at Glenn (formerly Lewis) 
Research Center in 1980s, culminating in an extensive compila- 
tion by Strazisar et al [25]. Shown in Fig. 6 is the test model 
of the rotor with 22 blades assembled. The rotor was designed 
for axial inflow and did not require inlet guide vanes, nor a stator 
stage. The design total pressure ratio is 1 .63 at a mass flow rate of 
33.25 kg/s (choked at 34.96 kg/s) and a rotating speed of 16,043 
rpm. Other geometrical dimensions and operating conditions can 
be found in the cited reference. The laser anemometry measure- 
ments acquired on streamsurfaces, starting at roughly one chord 
length upstream of the rotor and continuing through it till some 
distance into the wake, providing detailed data in the form of rel- 
ative Mach number and relative flow angle. Flow variables were 
also available at an upstream and a downstream planes. These 
data will be used to validate our computed results first, before 
building up the reduced order. 

The characteristics boundary conditions are employed for 
the inviscid boundaries: at the subsonic inflow boundary the left 
running (negative) Riemann variable is extrapolated from the in- 
terior domain; at the subsonic outflow boundary the static pres- 
sure is specified at the hub and radial pressure equilibrium as- 
sumed to calculate other radial points. The no-slip conditions 
are applied at the hub. But the shroud in this study is assumed to 
be an inviscid wall (or a streamsurface), this may cause some dis- 
crepancies of our results from the data, as will be remarked when 
appropriate. The mesh used for the CFD solution is shown for 
the blade tip in Fig. 7 in an overall view and two enlarged views 
showing the dense mesh at the blade surface and in the wake. 
Similar mesh distribution is also generated for all spanwise sec- 
tions. Two H-type mesh systems, one coarse and another fine, are 
used first to establish whether the coarse mesh is sufficiently ac- 
curate to be used for further aeroelastic analysis. The two meshes 
respectively consist of 77x43x45 and 104x63x80 grid points (re- 
sulting in 140,448 and 504,494 cells)-the three numbers respec- 
tively refer to the streamwise, blade-to-blade, and spanwise di- 



FIGURE 6: BLADE SHAPE OF THE NASA ROTOR 67 
AS MOUNTED ON THE HUB, TOGETHER WITH THE 
ALIGNED COMPUTATIONAL SURFACE. 


rections. This mesh density may be considered coarse in today’s 
CFD practice, however, Fig. 8 shows that the computed profiles 
of static pressure and total pressure and temperature ate the exit 
plane from both grids are essentially indistinguishable. These 
solutions are also comparable with the measured data and other 
published CFD results using finer meshes, for example in [5,22]. 
Hence we consider the coarse grid to adequate to provide suf- 
ficiently accurate aerodynamic forces to the structural analysis 
and thus, to be employed in this study. Moreover, our emphasis 
in this paper is to show the efficacy and validity of the proposed 
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(a) blade tip (b) leading edge 



(c) trailing edge 



(a) static pressure ratio 


(b) total pressure ratio (c) total temperature ratio 


FIGURE 8: PROFILES OF STATIC PRESSURE, TOTAL 

PRESSURE AND TEMPERATURE RATIOS AT AN EXIT LO- 
CATION WHEN THE ROTOR IS NEAR PEAK EFFICIENCY. 


model order reduction method for AE analysis. 

It is noted, however, that an overestimation is found in the 
static pressure ratio by the computation. This is probably caused 
by several simplifications committed in our computational setup: 
(1) we did not assume a boundary layer profile at the inflow 
boundary while in the experimental setup a solid surface is con- 
nected to the hub surface of the rotor, (2) the tip clearance is not 
taken into account and instead an inviscid slip wall is assumed at 
the casing, and (3) the hub wall is assumed adiabatic, hence pos- 
sibly giving rise to a higher temperature or pressure in the layer 
at the hub. This low-momentum layer at the inlet will continue 
to develop, growing through the rotor, resulting in a thickened 
boundary layer profile, in comparison with the computed result 
which indicates a fuller profile in a thinner layer. 




(c) 30% span from shroud 




(b) 10% span from shroud 




(d) 30% span from shroud 



(f) 70% span from shroud 


FIGURE 9: RELATIVE MACH NUMBER CONTOURS AT 
THREE SPANWISE SECTIONS, RESPECTIVELY 10%, 30% 
AND 70% MEASURED FROM SHROUD. 


The relative Mach contours at three span wise sections, re- 
spectively 10%, 30% and 70% measured from the tip, are com- 
pared for the peak efficient condition in Fig. 9, revealing the 
nearly normal shock wave across the blade passage at the tip sec- 
tion, but subsonic or low supersonic near the root. 

Finally, we plot the rotor 67 performance by the CFD solu- 
tion in comparison with the measured values, as shown in Fig. 
10. At the peak efficiency point, the solution gives a mass flow 
rate of 33.68 kg/s, a total pressure ratio of 1.651, and an effi- 
ciency of 0.9178. The calculated results are in good agreement 
with the data over the entire operating conditions. 
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TABLE 1 : MATERIAL PROPERTIES. 


Material 

Young’s modulus 

Poisson’s ratio 

Density 


(Pa) 


(kg/m 3 ) 

Titanium Alloy 

1.172xlO n 

0.3 

4539.5 


TABLE 2: MODAL FREQUENCY (HZ). 



1st mode 

2nd mode 

3rd mode 

Present 

369.8 

1009.4 

1622.9 

Doi [5] 

401.9 

1096.0 
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FIGURE 10: ROTOR 67 PERFORMANCE VS MASS FLOW 
RATIO: EFFICIENCY AND TOTAL PRESSURE RATIO. 


3.2 FLUTTER ANALYSIS FOR ROTOR 67 

The flutter characteristics depends on the structural proper- 
ties, in addition to the aerodynamic conditions. The material cho- 
sen for consideration is titanium alloy whose properties are given 
in Table 1, same as those used in [5,22]. (The material in Doi’s 
work was altered to give the Young’s modulus of 1 .422 x 10 1 1 (Pa) 
to place the first natural frequency away from the rotating fre- 
quency or its double.) 

The first three modal frequencies, calculated with the com- 
mercial software MSC/Nastran [12]onal5xl5 mesh, are listed 
in Table 2, along with the results by Doi. Despite using different 
values of Young’s modulus, these two predicted modal frequen- 
cies are quite close. These three mode shapes are shown in Fig. 
11; the first mode representing the bending, second mode the 
second bending, and third mode the torsion. 


We now apply the model order reduction technique es- 
tablished above to rotor 67 at the peak efficiency condition. 
The time-dependent aerodynamic force is built using the above 
Volterra series with the first mode displacement given as: 


= 5.0x10 5 sino)f (23) 

where (O is the first natural frequency of the structure. 

It is noted that the small amplitude is chosen in Eq. (23) to 
ensure linearity assumed for the current ROM formulation. Since 
the flutter boundary estimated by the linear theory is independent 
of the perturbation magnitude, it is not critical what value is used 
as long as the value is small. The time step used in the time 
integration is chosen to be sufficiently small that time accuracy is 
maintained; in this study the time step is 2.0 x 10 — 5 (^), allowing 
about 30 time-intervals in the period of the highest frequency 
mode considered. 

In Fig. 12, we validate the accuracy of the ROM-CSM 
model (system of Eqs. (17) and (7)), by comparing the first 
three modal forces of rotor 67 blade. The aerodynamic force 
is obtained either by solving the full Reynolds-averaged Navier- 
Stokes equations or by using the ROM with 150 degrees of free- 
dom (i.e. the order of matrix A a in the state space model) for 
an input defined by Eq. (23). It reveals that the third mode 
(torsion) is most dominant and the weakest is the second mode 
(second bending). The close agreement between the full and 
reduced models confirms the accuracy of the current Volterra- 
series-based ROM. Discrepancy is seen in the second mode, but 
this mode is less important than the other two. 

The aeroelastic ROM system consists of 150 degrees of free- 
dom in the aerodynamic ROM and 6 in the structural represen- 
tation (displacement and velocity), thus resulting in 156 DOF in 
total for the entire AE ROM. The instability critical point can 
be determined by increasing total density at inlet boundary. As 
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(a) first mode: bending 


(b) second mode: second bending 


1 st mode 



0.00 0.01 0.02 
time(s) 


FIGURE 12: COMPARISON OF MODAL FORCE OBTAINED 
BY THE FULL ORDER AND REDUCED ORDER SOLU- 
TIONS. 



(c) third mode: torsion 


FIGURE 11: MODE SHAPES OF THE ROTOR 67 STRUC- 
TURE AND INTERPOLATION BETWEEN THE STRUC- 
TURAL AND AERODYNAMIC GRIDS. BLACK: NON- 
DEFORMED GRID; BLUE: STRUCTURAL GRID; RED: 
AERODYNAMIC GRID. 



Re 


FIGURE 13: EIGENVALUES OF THE 156-ROM. 


shown in Fig. 13, the 1st structural mode eigenvalue crosses the 
imaginary axis, i.e., the eigenvalue becomes positive, indicating 
an amplification of structural displacement. Figure 14 displays 
the blade displacement predicted by the 156 aeroelastic ROM at 
the flutter condition, the third and first modes are the two most 
dominant ones while the second mode is nearly negligible. The 
dynamic pressure needed to induce flutter is qoo = 1.455vl0 6 Pa, 
nearly 10 times larger than the baseline operating condition at 
= 1.416vl0 5 Pa. Hence the rotor made with the material spec- 


ified in Table 1 is determined to be structurally stable under the 
chosen operating condition, with a high margin of safety, when 
only an isolated blade is considered, this finding consistent with 
that in [5,22]. However, blade row interactions, such as the effect 
of upstream inlet guid vane, can induce forced vibration in rotor 
blade, thus altering its flutter characteristics, see study in [19]. 
The ROM strategy presented here can also serve as an efficient 
and reliable way of investigating the effect of inter-blade interac- 
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displacement(m) 


tions. 

Finally we remark on the primary motivation of employing 
ROM, while under the foremost requirement of preserving ac- 
curacy. For performing an aeroelastic analysis over a complete 
sinusoidal cycle (Eq. (23)), the full-order (CFD-CSD) model 
takes 10.8 hours on a Xeon(R) W3530 computer with Intel(R) 
Compiler compared to 0.56 seconds used by the ROM, a whop- 
ping savings by over 19,200 times. This shows the tremendous 
value of using the ROM when searching for the flutter bound- 
ary shown in Fig. 14, or when conducting design optimization, 
both of which will otherwise require enormous computational 
resources. 


1 st mode 



FIGURE 14: FLUTTER RESPONSE IN TERMS OF DIS- 
PLACEMENT OF THREE MODES WHERE THE TORSON 
AND BENDING MODES ARE DOMINANT AND THE SEC- 
OND BENDING MODE IS MINIMAL. 


CONCLUDING REMARKS 

We have presented an accurate and efficient method for per- 
forming aeroelastic analysis of a modern transonic compressor 
blade, NASA rotor 67. The CFD code, using k-co-SST tur- 
bulence model and AUSM + -up numerical fluxes, for providing 
aerodynamic forces has been validated against measured data. 
The structural motion based on finite element analysis is coupled 
with fluid motion, the coupling is further modeled by the state- 
space representation to achieve considerable reduction in compu- 
tational cost, while preserving the solution accuracy. The linear 
state-space system is formulated by keeping only the first-order 
Volterra kernel. The obtained reduced order model is shown to 
be in excellent agreement with the full (original) model. Hence it 
can be employed to provide an effective aeroelastic analysis tool, 
specifically for defining the flutter boundary. 
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